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The r~ n interaction energy, 1 < n < 3, for a infinitely periodic system with 
explicit charges and a neutralizing, uniform background charge density is derived. 
An Ewald based expression for this energy has an extra term proportional to the 
square of the total explicit charges of the system. This expression may be useful for 
simulations in which explicit charge neutrality does not hold or for which the total 
explicit charges is a fluctuating quantity. 
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1. INTRODUCTION 

Simulation of atomic systems very often entails the calculation of the Coulomb inter- 
action of point charges. The commonly used periodic boundary conditions embed the 
simulation box in an infinitely periodic system. The interaction energy per unit cell volume 
of this infinitely periodic system, of N charges Zi,i = l, . . . , N, can be written as 

N N 

" = EEE' i r T ii" * (1) 

1 i= i j= i i i r * r J i \ 

Here 1 represents the lattice vectors defined by the simulation box, and the prime on the 
sum excludes the term 1 = when i = j. The Coulomb interaction corresponds to the case 
n = 1. To make the discussion and derivation slightly more general, the inverse power is 
allowed to vary in the range 1 < n < 3. It has been shown that when 1 < n < 3, the 
sum converges only when the overall neutrality is met, i.e., J2i=i z i = UL F° r n = 1, 
the sum is conditionally convergent, i.e., the value of the sum depends on the order in 
which the terms are added. An elegant, conventional way to calculate this sum, the Ewald 
summation, introduces a convergence function that separates the sum into two components. 
The first component converges rapidly in the real space and the second converges rapidly 
in the Fourier space [2, 3, 4, 5]. The k = term, k being the reciprocal lattice vectors, in 
the Fourier space depends on the shape of the infinitely periodic system [6, 7] and can be 
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recast as a surface integral [8]. Its value can be determined analytically in simple cases or 
evaluated numerically in the general cases. Ignoring this term corresponds to the boundary 
conditions with an infinite dielectric constant [5]. Faster algorithms for the calculation 
of the Coulomb interaction have been developed, and references can be found in a recent 
review paper [9]. 

Problems arise during simulation of charged systems when some ions are treated im- 
plicitly. For example, when a peptide with basic or acidic amino acids is embedded in a 
shell of water molecules, the explicit charges may not meet the neutrality condition and, 
hence, result in an infinite energy. One possible way to overcome this problem is to add 
counter ions or fictitious charges in a way to maintain charge neutrality [10]. This paper 
presents an alternative way to arrive at a finite energy, by adding to the system a uniform, 
neutralizing charge density p = — Z i/V c , V c being the volume of the simulation 
box, or unit cell. This background charge density adds negligible amount of computational 
time, and the result obtained from this derivation may be useful in practical simulations. It 
allows for a computationally convenient approach to ionic systems where some ions and 
atoms are treated implicitly. In particular, it facilitates simulations of ionic systems in the 
grand canonical ensemble, where charged molecules can be added to or removed from the 
system. The accommodation of the inverse power n = 2 may provide useful expressions 
for screened Coulomb interactions. For example, one way to model the solvent effects in 
the AMBER forcefield [11, 12] for biological systems is use a distance-dependent dielec- 
tric constant such as e/eo = 4r, where r is in A [13]. Therefore, the screened Coulomb 
interaction is proportional to r~ 2 and can be calculated by the results developed for n = 2. 

This paper is organized as follows. In Sec. 2 a Ewald based, computationally useful 
expression for the energy of this system is derived. It is shown that that an extra term 
proportional to {J2iLi z i) 2 ar i ses when the explicit charge neutrality does not hold. The 
result and the k = term originating from the boundary conditions are discussed in Sec. 3. 
The conclusion is given in Sec. 4. 

2. DEVELOPMENT 

Consider in a infinitely periodic system the potential tpi contributed by a point charge 
Zi located at r i7 and its neutralizing uniform density pi — —Zi/V c . For 1 < n < 3, the 
potential can be written as 



V>i(r) = lim 

L^oc 



V 1) - — / -, r ${L- v')dr' 

^ |r-ri-l|" 1 ; V C J |r-r -r'|" 1 ' 



(2) 

where L is the characteristic dimension of the periodic system and 

<KL- I) = I 1, re V{L) (3) 
1 ' ' \ 0, r £ V(L) W 

is a cutoff function that is added to specify the shape of the periodic system. The vector r 
represents the geometrical center of the cell, r = J v rdr/V c . The procedure of including 
the shape of the volume followed by taking the limit is non-trivial and leads to the k = 
term in the Coulomb interaction case [8]. Note that when n > 3, tpi(r) diverges due to 
the singularity at r' = r — r . Next, we introduce a convergence function f(r), which 
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tends rapidly to zero at infinity. This convergence function is used to split ^(r) into two 
components: 



^! K) (r) = lim 

L — >oo 



^ |r- ri -l| n V ; 

14 7 |r-r -r'|" v ' y 
l-^dr-Fi-ll) 



ip\ (r) = lim 

L^oo 



E 



Ir-r, : -l| 



0(L; 1) 



r - r - r' 



(4) 



Here (r) and ^- ; (r) stand for the short-ranged and long -ranged parts of the sum, 
respectively. Assuming the convergence function <p(r) decays fast enough so that each 



term in r/;> J (r) is finite, the limit i — > oo can be formally taken and 

^ ( j " ^ |r- ri -l|» Kei |r-r -r''" " : 
<£>(|r -r; - 1|) /• p(|r'|) , 



E 



Ir-ri-ll 



(5) 



Note that the average of ipl over the simulation box is zero: 

V* (r)dr = k z ^T^TF dr -J v vJ Tr drdr 

J 1 k'l 



¥>(|r'|) , , 

-ctr — Zi 



0, 



(6) 



where J y denotes integration over the simulation box. The long -ranged J (r)isbest 
evaluated by rewriting it as a sum of Fourier components. Let [/(r)] denote the Fourier 
transform 



F k [/(r)] = y"/(r)exp(-ikT)dr. 



(7) 



The Fourier components of ^ F ^(r) are non-zero only at the reciprocal lattice vectors k. 
Therefore, 



l-^(|r'|) 



exp [ik ■ (r -Ti)} + Wj(r), 



(8) 



where 



Wi(r) = lim 



l-^(|r-ri-r'|) 1 - <p(|r - r - r'| 



r - r, : - r' 



r - r - V 



0(L; r')rfr' 



(9) 
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For n = 1, u>i(r) leads to the non-vanishing k = term as L — > oo. However, it 
diminishes when n > 1. We shall neglect this term for now and discuss it later in Sec. 3. 
Uses of Eqs. (3), (4), and (8) give 



V>;(r) 



E 



(p(\r 



vX F > 



r-r;-l|" K 7 |r'|™ 
l-^(|r'|))- 



exp [ik • (r - r-j)] . 



(10) 



Now we add up the contributions from all charges and the background density and let the 
total potential be ip(r) = J2iLi *Pi( r )- The total interaction energy, per unit cell volume, 
is therefore 



U 



1 N 



V'(r) 



Zi 



-p / ip{v)dv. 



(11) 



The bracketed term in Eq. (11) stands for the potential acting on the point charge i, and 
can be rewritten as 



V>(r) 



r - r, 



JV 

EE- 

i=i l 



'z M\r-j -tj -lj 
ki-r, -II" 



l-^(|r|) Ef=i^ 



-Z< lim ,7' [/ - J ^hr^dv' 

IrHo r" y c 



7 |r'|" 



4l> 



i-v(M) 



exp(ik-ri)5(k) (12) 



Here S*(k) = EjLi z j exp(— zk • rj) is the structure factor. If we ignore the surface 
contribution, the k = term is zero, as is the second term at the right hand side of Eq. (1 1). 
Finally, substituting Eq. (11) into Eq. (11) yields an expression for the total energy: 



JV N 



2 i=l 3=1 1 

T N z- 



2Vr. 



lim 

2 |r|-o r 



i-v(M) , i 



7^ 

7 |r'| 



+ 



2V C 



£|S(k)| 2 F k 



i-v(M 



(13) 



The convergence function <p(r) should be chosen in a way that each term of Eq. (12) is 
finite and computable. One possible choice is tp(r) = I ^r(n'/2)' ^ where F(n/2) and 
r(n/2, ar 2 ) are the gamma function and the incomplete gamma function, respectively. 
The parameter a adjusts the contributions from the real and the reciprocal parts and should 
be optimized for the system of interest. With this choice of ip(r), the Fourier transform of 
(1 — <^(|r|))/|r|™ has been shown to be [3] 
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and 



l-<p(\r\) 2a"/ 2 



Km 
M-o |r| 



nT{n/2)' 



(15) 



The short-ranged part of the contribution from the background charge density can be 
evaluated analytically to be 



¥>(|r|) 



dr 



9^3/2 



|r|" (3-n)r(n/2)a( 3 -")/ 2 ' 
Therefore, the interaction energy can be written as 

N N 



< n < 3. 



(16) 



U = 



1 



2T(n/2) 



EEE" 2 ' 1 '^ a|ri - rj- - 1| 2 ) 



i=i j=i l 



lr ■ — r — 1 1™ 



N 



2a™/ 2 



N 



+ 



(3 - n)y c a( 3 -™)/ 2 ^ Zl J n E 

2 3 -"tt 3 / 2 ^ |S(k)| 2 /3-n fc^\ 



(17) 



In the Coulomb interaction case, n = 1, cp(r) is the complementary error function 
erfc(a 1 / 2 r) and Eq. (16) is reduced to the Ewald sum [5] plus an extra term proportional 

to (Ef =1 ^f-- 

tt 1 V , V^Y^' erfc(a 1 / 2 |r i - Tj - 1|) tt /A \ 

^=1 = 2EEE lr-.--v.--ll 9V^(E*J 

»=i j=i l 



|r* — rj- — 1| 



(-) E^ + ^E 



2tt ^ |S(k)| 2 / fc 2 \ 



k 2 



CX P - V • 



i=l k^O 

For n = 2, ^>(r) = exp(— ar 2 ) and the corresponding equation is 



4a / 



(18) 



AT N 



u n= 



^EEE' 

i=i i=i i 



exp(-a|r; - r,- - 1| 2 ) 7r 



3/2 



AT 



Ti-Tj -1| 



1/2 E ; 



K ^ k 

k^O 



2aV2 



(19) 



3. DISCUSSION 

The k = term described by Eq. (8) does not diminish when n < 1. Let w(r) = 
Si=i w i( r )- When w(r) is multipole-expanded and the limit L — > oo is taken, the first 
term becomes zero due to charge neutrality. Thus, the leading term in w(r) is the second 
term, which is non-zero only at the surface of the periodic system and can be recast as 

f\ v n f (r - r M ) • r'd • n 
w(r) - lim — / — 2 dS, (20) 

l^oo V c J S(L) r' 



6 



! Please write \authorrunninghead{<Author Name(s)>} in file ! 



where 

Sill N r » + I HiLl Z i\ V nn 

Ei=i N + IE i= i^l 

is the "center of mass" of the simulation box, d = EiLi z i( r i — r o) is the dipole moment 
of the simulation box, and S(L) is the surface of the volume with a characteristic dimension 
L. Note that both the explicit charges and the background charge density contribute to the 
dipole moment. It is interesting to see that u>(r) is the only term that contributed to the 
average potential over the simulation box, i.e., 



n 



(r - r M ) • r'd • n 



I™ T7 — ^To. dS. (22) 



Vc J S{ L) r'"+ 2 

Note that the average of both the short-ranged part and the k ^ terms over the simulation 
box is zero. 

The contribution of w(r) to the energy per unit cell is 

N 

(23) 



1 % ■> 1 f 

U s = ~Y^ z M r i) - -p uj(r)dr. 



Substituting the expression for u(r) into Eq. (23) gives 



n 



d • r' d • n 



U B = lim — / dS. (24) 

l-,oo 2V C J S(L) r' n+2 

Clearly, this term diverges when n < 1 and diminishes when n > 1 as L — > oo. For n = 1 
it contributes a finite energy to the system [8]. It is interpreted as the interaction between 
the dipole moment of the system and the surface charges at the boundary. Ignoring this 
term corresponds to using boundary conditions with an infinite dielectric constant. 

It is worth mentioning that, for the n > 3 case, an expression having the same form as 
Eq. (16) was derived in ref. [1] without the addition of a neutralizing background charge 
density. Here the same equation is extended to the case 1 < n < 3 with a different physical 
interpretation. 



4. CONCLUSION 

The r~ n interaction energy, 1 < n < 3, for a infinitely periodic system with charges 
and a neutralizing, uniform background charge density is derived. Following the Ewald 
transformation, the interaction energy is split into a real space sum and Fourier space sum, 
and both are rapidly convergent. The Ewald based expression for this energy has an extra 
term proportional to the square of the total explicit charge of the system. This expression 
may be useful for simulations in which explicit charge neutrality does not hold or for which 
the total explicit charges is a fluctuating quantity. 
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